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ABSTRACT 


Parameter estimation methods were used to obtain estimates of stability and control derivatives 
for the Marine Corps BQM-147 Unmanned Air Vehicle. The results from a simple, PC-based linear 
model and those from a more robust non-linear model, pGst, were compared. A Cramer-Rao bound 
was used to assess the accuracy of the estimates for both methods. The bounds were high for both 
the longitudinal case and the lateral-directional case due to the limited maneuvers tested, high levels 
of noise in the same general frequency range as the control input, and the lack of body-angle data. 
The linear model failed to provide estimates for the lateral-directional case. Though the results may 
be used as starting points for a dynamic model of the aircraft it is recommended that the flight test 
procedures be modified to address the issues raised concerning noise, recorded signals and the need 
for repeated maneuvers. 
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I. INTRODUCTION 

Symmetrical-wing, delta-planform (SYMDEL) unmanned air 
vehicles (UAVs) have been in development and have flown in a 
variety of mission profiles since the early 1970's. This 
delta-planform, symmetrical-airfoil configuration was designed 
to provide equal performance in range, speed and stability in 
either an upright or inverted attitude [Ref. l:p. 23]. With 
this attribute the vehicle can readily accommodate a variety 
of payloads which makes it ideally suited for specific 
military applications or as a testbed vehicle for new 
component technologies. 

The EXDRONE UAV was developed in the early 1980's at the 
Johns Hopkins University, Applied Physics Laboratory (APL). 
The vehicle design is a descendent of the earlier SYMDEL 4 
with a constant-chord section inserted into the middle of the 
delta planform to provide an accessible avionics bay [Ref. 
l:p. 25]. The original APL version of the EXDRONE, referred 
to as the baseline configuration throughout this document, was 
intended to meet a United States Marine Corps (USMC) 
requirement for a UAV to carry an expendable VHF 
communications jammer. Takeoff, climb out and turn to the 
desired heading are accomplished by radio control. Once the 
mission heading is achieved, a radio control command switches 
control of the vehicle to the autopilot. The autopilot then 
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controls the vehicle for the duration of the mission. The 
name EXDRONE implies that the vehicle was intended to be 
operated autonomously and to be expendable at the end of its 
mission; however, radio control can be reacquired if necessary 
and the vehicle can be recovered. 

In preliminary flight tests the baseline configuration of 
the EXDRONE was found to be difficult to handle for less 
experienced pilots. The vehicle had weak stability and 
control characteristics in low-speed flight and had a tendency 
to depart near stall [Ref. 2:p. 1]. The primary mission of 
the EXDRONE is to dash outbound at a speed of about 100 miles 
per hour to a predesignated target area and then to loiter in 
the area and jam enemy VHP communications. With low speed 
flight being an equally important part of the EXDRONE mission, 
the USMC UAV Project Office contracted with the NASA Langley 
Research Center, Flight Dynamics Branch, to conduct 
exploratory wind tunnel and flight tests in order to determine 
the aerodynamic characteristics of the vehicle and to provide 
modifications for aerodynamic improvements [Ref. 2:p. 1]. As 
a result of the wind tunnel study several modifications were 
made to the baseline configuration in order to increase pitch 
control for nose-up trim, improve the stall departure 
characteristics, minimize the adverse yaw characteristic and 
reduce the minimum level flight speed as a means of increasing 
endurance. It was this modified configuration EXDRONE that 
became the subject of the flight test portion of the NASA 
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Langley study. A detailed description of the wind tunnel 
tests on both the baseline and modified configurations can be 
found in Reference 2. 

The modifications to the wing tips and control surfaces 
that are described in Reference 2 were included in the 
production contract for the BQM-147A. The following 
additional modifications are also a part of the production 
configuration. 

1) Position and orientation of the engine 

2) Shape of the nose/cowl 

3) Number and type of skids 

In October of 1991, Developmental Testing (DT) was 
conducted on the BQM-147A at the Dugway Proving Grounds, Utah. 
The purpose of the tests was to determine the suitability of 
the vehicle in the VHF communications jamming mission and to 
determine the effects of the communications jammer payload on 
the handling qualities of the airframe. The BQM-147A 
control surfaces are manipulated by electrical servos that are 
controlled by an autopilot. The autopilot receives command 
signals from either a modular radio control unit for direct 
pilot input or from an onboard microprocessor containing a 
pre-programmed mission scenario for the autonomous mode. The 
autopilot processes the error signal from a rate gyro to 
augment roll damping in flight [Ref. 3:p. 20]. During DT the 
autopilot had to be fine tuned by hand in a hit or miss 
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fashion because the vehicle control laws had not been well 
defined during the pre-production work. It was at this point 
that the USMC UAV Project Office decided to fund a series of 
flight tests from which time history data could be used to 
determine the aerodynamic characteristics of the aircraft. 

Research being done at the Naval Postgraduate School 
(NPS), Department of Aeronautics and Astronautics, using UAV 
testbeds has been well received by the UAV Joint Project 
Office. Specifically, aerodynamic parameters obtained by 
numerical methods and by wind tunnel tests [Refs 4 and 5] have 
been used by the Target Simulation Lab, Naval Air Warfare 
Center, Weapons Division (formerly Pacific Missile Test 
Center) , Pt. Mugu, CA., in the development of a math model for 
the Pioneer UAV flight training simulator. More recent thesis 
work by U.S. Navy LCDR Robert Graham incorporated a personal 
computer (PC) based parameter estimation capability into the 
NPS research program [Ref. 6]. Being advised of this desktop 
parameter estimation capability, the USMC UAV Project Office 
expressed an interest in an application of this routine using 
measured data obtained from the next series of flight tests. 

The objectives of the application listed in Reference 7 are: 

1) To document the performance of the aircraft. 

2) To reprogram the autopilot on the production aircraft. 

3) To form a basis for suggesting changes to the production 
aircraft which will improve stability and performance. 

4) To develop a PC-based flight simulator which will be 
used to train operators. 
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Aircraft parameter estimation routines are used to 
estimate aircraft stability and control derivatives from 
actual time history data. The stability and control 
parameters extracted from these routines provide an accurate 
model of the actual vehicle and can be used in analyzing the 
effects of proposed control system changes on aircraft 
stability and handling qualities. Parameters obtained from 
observed data can also be used to provide an accurate model of 
the vehicle for greater fidelity in flight simulators. 

The purpose of this study was to use the flight test time 
history data found in Appendix B in a parameter estimation 
model to estimate both the longitudinal and lateral- 
directional stability and control derivatives. Two parameter 
estimation models were used. The first was an application of 
a PC-based parameter estimation program by LCDR Robert Graham 
[Ref. 6]. This model was selected because it is readily 
available, it operates in a familiar DOS environment and it 
generally provides good results when compared with more 
complex mainframe programs [Ref. 8J. The disadvantage of this 
first parameter estimation routine is that it is a linear, 
maximum likelihood estimator which assumes that the actual 
system satisfies the assumptions that were made in the 
linearization of the model. The second model, named pEst, is 
an interactive program for parameter estimation in non-linear 
dynamic systems. Incorporation of the pEst routine into the 
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ongoing flight test research at NPS provides a capability to 
model the dynamic behavior of aircraft that cannot be 
appropriately modeled using the assumptions inherent in the 
shorter, PC-based MMLE-3 code. Because of suspected non- 
linearities in the more extreme areas of the flight envelope 
of the EXDRONE, the linear model may not be sufficient and the 
pEst routine might better capture the behavior of the system. 
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II. BACKGROUND 


A. HISTORICAL REVIEW 

The concept of using remotely piloted vehicles (RPV's) and 
unmanned aerial vehicles (UAV's) in a tactical combat role had 
limited acceptance within the military community from the time 
of its first proposal in 1915 until early in the last decade. 
The RPV is controlled by a human operator by means of a radio¬ 
control guidance link using either direct operator visual 
contact or onboard video imagery to sense flight parameters. 
Initial RPV testing was limited by pre-WW II radio-control 
technology but improvements in guidance technology led to 
successful developmental flight tests of ordnance carrying 
RPV's prior to the attack on Pearl Harbor in December 1941. 
There was limited use of these lethal RPV's throughout the 
latter stages of World War II and the Korean Conflict. 
Widespread use of these platforms and full-scale development 
of follow-on vehicles with greater combat potential were met 
with opposition by military planners who preferred to use 
scarce resources on more expensive and ostensibly more 
technically-challenging guided-raissile programs [Ref. 10]. 

Advancements in radio-control and digital communication 
technologies coupled with an operational need to defeat more 
modern military anti-air warfare weapons led to an expanded 
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role for a new generation of non-lethal UAV's. The UAV 
differs from the RPV in that the guidance signal, from either 
a human operator via radio-control uplink or from on-board 
navigation subsystems, is input to an auto-pilot which 
generates the appropriate control deflection. UAV's are 
generally known for their ability to operate autonomously over 
the entire range of a preprogrammed mission profile. During 
the Vietnam conflict the Ryan Model 147 Firebee was used for 
photo reconnaissance, electronic intelligence gathering, 
communications intelligence and covert psychological warfare 
missions [Ref* ll:p. 20], In spite of the tremendous success 
of UAV's in Southeast Asia, continued support for further UAV 
development waned in the wake of higher priority military and 
domestic programs. Not until the unprecedented and highly 
successful use of UAV's by the Israelis in the Bekaa Valley in 
1982 did U.S. military planners become convinced of the 
tactical usefulness on non-lethal UAV's [Ref. 12]. 

Recognizing the need for common and interoperable systems, 
Congress in 1988 directed the Department of Defense (DoD) to 
consolidate the management of non-lethal UAV programs. The 
DoD established the UAV Joint Project Office (UAV JPO) and 
appointed the United States Navy as the executive service. 
The mission of the Joint Program Office is to expeditiously 
field quality UAV systems which provide a significant tactical 
advantage to operational commanders. The objective of the UAV 
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JPO is to field a family of interoperable systems that 
optimize commonality [Ref. 13:p. 7]. 

During Operation Desert Storm, UAV's provided tactical 
commanders an operational capability that was not previously 
available. The Marines operated three types of UAV's in 
various mission roles that included real time reconnaissance, 
surveillance, target acquisition, spotting and bomb damage 
assessment. When operation of the Pioneer UAV was jeopardized 
by parts shortages and while the more fragile Pointer was 
unable to fly in winds greater than 15 knots, the Marine Corps 
augmented their UAV operations by reassigning a number of 
EXDRONE's from a research and development program to be 
operationally deployed from airfields near the Saudi Arabia- 
Kuwait border [Ref. 14]. These vehicles were reequipped with 
color video cameras and were used to find attack routes 
through enemy defenses. 

This demonstrated capability of UAV's to perform multiple 
functions and the ability to provide near real-time 
intelligence in a high-threat environment without loss of life 
has affirmed the commitment by the Services to integrate UAV's 
into the force structure. The 1992 DoD Unmanned Aerial 
Vehicles Master Plan [Ref. 13] indicates that the planned 
procurement of UAV's is in excess of 2300 vehicles. The core 
strategy of the UAV JPO is the establishment of a family of 
UAV systems which indicates that these systems must allow for 
growth in performance and must readily accommodate changing 
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payloads. As a part of this core strategy the JPO has 
initiated a Very Low Cost (VLC) UAV program. These VLC 
systems are used to demonstrate and evaluate the utility of 
UAV's in tactical units. 

In November 1991, the UAV JPO awarded a contract to BAI 
Aerosystems, Inc. for the production of 110 EXDRONES. At the 
present time, operational test and production for the Marine 
Corps Jammer program is not funded and the EXDRONE vehicles 
are being used in the JPO VLC program [Ref. 13:p. 49], If the 
production option in the BAI contract is exercised, it is 
expected that as a result of the VLC program a new effort to 
improve range and payload capacity will begin. Having 
accurate estimates of stability and control derivatives will 
allow designers to incorporate good flying qualities from the 
earliest stages of design and avoid repeating errors from 
earlier versions of the model. 


B. PARAMETER ESTIMATION 

In Reference 15 Ljung points out that the construction of 
a model from observed data involves three basic entities: 

1) The data record. 

2) The set of candidate models. 

3) A rule by which candidate models can be assessed using 
the data. 
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In an application to aircraft stability and control, the 
flight test maneuvers should be designed so that the data 
collected will provide the maximum value to the analyst. 
Consequently it essential that the flight test program have a 
clear set of objectives so that the flight test maneuver will 
cover the area of the envelope that is to be explored. 
Instrumentation must be considered so that the appropriate 
parameters are measured to a degree of accuracy that will be 
suitable to the objectives of the program. The set of 
candidate models must also be considered. The equations of 
motion derived from analysis of forces and moments acting on 
a typical aircraft are well defined. For most applications a 
linearized set of the equations is adequate for analyzing 
aircraft stability and control. If a linear model is all that 
is available for the analysis, then only small- perturbation 
maneuvers should be included in the flight test. If envelope 
expansion is the objective of the flight test, then maneuvers 
that gradually expand the envelope may provide data that 
indicate a non-linear model is necessary. The rule by which 
the models are assessed is generally determined by how well 
the model is able to reproduce the measured data. For 
example, the modified maximum likelihood program uses an 
output-error approach that minimizes the difference between 
the measured parameters and the model estimate of the 
parameters in successive iterations. This construction of a 
mathematical model of a dynamic system from observed data is 
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know as System Identification [Ref. 15:p. 1]. 

Estimation of aircraft stability and control derivatives 
from flight test time history data is one of many practical 
applications of the basic system identification techniques. 
Of particular note in the aeronautics industry is the work of 
Lawrence Taylor, Kenneth Iliff, Richard Maine and James Murray 
of the NASA Ames Research Center, Dryden Flight Research 
Facility. Their approach to the estimation of aircraft 

stability and control derivatives through the application of 
system identification techniques can be found in References 9, 
17 and 18. The method of application can be summerized in the 
following steps: 

1) Formulate the dynamic equations that describe the model. 

2) Identify from the form of the dimensional derivatives 
the parameters whose values are unknown. 

3) Collect flight test data. 

4) Infer the values of the unknown parameters by adjusting 
the values iteratively until the computed response best 
matches the measured response. 

Throughout the 1980's the most widely accepted aircraft 
parameter estimation routine was the Modified Maximum 
Likelihood, version 3, written by Iliff and Maine. The MMLE-3 
program evolved from an earlier Newton-Raphson parameter 
estimation routine by Taylor and Iliff when a more versatile 
program capable of handling larger amounts of data was 
required by the industry. When applying system identification 


12 





techniques to a real world system, modeling errors are 
unavoidable since it is likely that the real system does not 
completely satisfy all of the assumptions on which the model 
was developed. The existence of unmeasured inputs along with 
these modeling imperfections perturb the states and to at 
least some degree invalidate the results of the identification 
algorithms [Ref. 15]. The advantage of MMLE-3 is that the 
algorithm is capable of applying the system identification 
techniques in the presence of state noise and thus produces a 
more optimistic error bound estimate. Full details of the 
MMLE-3 program can be found in References 17 and 19. The 
minimization routines and the criterion used in a personal 
computer based adaptation of MMLE-3 are discussed in more 
detail in later sections. 

For many aircraft applications the linearized equations of 
motion provide a reasonably accurate model of the system for 
small disturbances. For some applications however, the linear 
model is not sufficient and a more accurate non-linear model 
is required. Flight testing at the Dryden Flight Research 
Facility highlighted the limitations of the linear model when 
applied to maneuvers at extreme flight conditions and with 
some unique aircraft configurations that exhibited non-linear 
dynamic behavior. In 1987 Murray and Maine introduced a 
program for parameter estimation in non-linear dynamic 
systems. The capabilities of the pEst program are described 
in Reference 9. 
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Historically the wind tunnel has served as the primary 
source of aeronautical data for flight vehicles. The 
advantage of the wind tunnel is that flow conditions can be 
accurately controlled, but other influences such as wall 
interference are difficult to account for. The wind tunnel is 
also limited in its ability to achieve dynamic similarity with 
actual flight test conditions. There are numerous reasons why 
wind tunnel data alone are not sufficient and often the 
results are validated by flight test. In Reference 18, Maine 
and Iliff point out that there is a tendency in practice to 
emphasize positive results and to down play the role of flight 
test in validating wind tunnel data or analytical predictions. 
A discussion of the errors that are found and corrected 
because of the flight test are often omitted from the 
published reports. Maine and Iliff go on to say that the 
utility of stability and control parameter estimation can best 
be evaluated by comparing the predicted data prior to the 
first flight test with the best estimates from combined flight 
data and predictions at the conclusion of the flight test 
program [Ref. 18:p. 4]. The task then is to take maximum 
advantage of previous information, both analytical and 
experimental, and to use the estimation routines to fill the 
gaps. The parameter estimation routines used in this study 
use a-priori knowledge of the system under investigation to 
provide to the user the most accurate estimates of aircraft 
stability and control derivatives. 
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III. 


DEVELOPMENT OF THE MATHEMATICAL MODEL 


The analysis of a dynamic system begins with the 
derivation of a mathematical model that is reasonably 
accurate. Typically the mathematical model is a set of 
differential equations obtained by the application of physical 
laws that govern the system. The response of the system to a 
given input may then be obtained by solving the differential 
equations. 

Development of the equations of motion begins with an 
application of Newton's second law to a rigid body in a non¬ 
rotating inertial-axis system. Motivation is given to make 
the transformation to a rotating, body-axis reference system. 
Euler angles are introduced and developed in some detail as a 
means of referencing body axes to the inertial axes. Finally 
a brief description of the aerodynamic model is presented in 
the form of non-dimensional force and moment equations. 

A. APPLICATION OF PHYSICAL LAVS 

Newton's second law states that the acceleration of an 
object is proportional to the net force exerted on the object. 
This law leads directly to the fundamental concepts of the 
conservation of linear momentum and the conservation of 
angular momentum [Ref. 20:p. 143]. In a non-rotating inertial 
reference frame the conservation laws can be expressed in 
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terms of the vector equations: 

uv) (i) 

(2) 

These equations state that the summation of vector forces (F) 
acting on a body is equal to the time rate of change of linear 
momentum (mv) and the summation of external moments (M) acting 
on a body is equal to the time rate of change of angular 
momentum (H). 

The externally-applied aerodynamic forces and moments 
acting on an airframe are due primarily to airflow conditions 
and control-surface deflections. The force and moment 
equations can be broken down into vector components along the 
longitudinal (x) , lateral (y) and vertical axis ( 2 ) of the 
body with respective velocity components u,v and w. 


d 

dt 

(mu) 

F y =-f c (mv) 


(3) 

L = 

-4- H x 

dt x 


S5 

11 

(4) 


Figure 1 illustrates the convention for the forces and moments 
in a body-fixed reference frame about the center of gravity. 
The components of the velocity vector V are u, v and w aligned 
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along the x, y and z axes respectively. The moments due to 
aerodynamic forces are L, M and N. 



References 21 and 22 detail the development of the moment 
equations referred to a moving center of mass. The result is 
a vector expression for angular momentum in component form 
given by: 
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/ 

piy 2 * z 2 ) 


g(xy) +r(xz) 


H y 

-f 


g(x 2 + z 2 ) 

- 

pixy) +r(yz) 



* vOi 

< 

r(x 2 +y 2 ) 


p[xz) + g(yz) . 

> 


dm 


(5) 


The expression for the angular momentum may then be written 
as: 

* = [*m] O (6) 

where the vector n contains the components of angular velocity 
p, q and r. The matrix [I m ] is a symmetric matrix of the 
moments and products of inertia which can be expressed by: 



t 

-I 

-I 


X 

xy 

xz 


~^xy 


-Ixz 


-*xz 

~!yz 



(7) 


with components: 


I x ~ J (y 2 + z 2 ) dm 

Iy~f (x 2 + z 2 ) dm 
J = f (x 2 +y 2 ) dm 


Ixy* fjxy) dm 
I xz = [ (xz) dm 

J V 

Z yz =fj y z) dm 


( 8 ) 


With the inertia matrix and angular velocity vector 
defined, equation (2) can be rewritten as: 

Equation (9) indicates that in a fixed reference frame, if the 
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airframe rotates then the matrix I ro will vary with time. By 
making the transformation to a rotating coordinate system that 
is fixed to the airframe, the inertia matrix becomes time 
invariant and the analysis becomes much simpler [Ref 21]. 

In making the transformation from the inertial reference 
to the rotating body reference, the vector equations for P and 
M can be expressed as: 


F = m — (v) x (mv) 
dt 


( 10 ) 


M = 


d 

dt 


(H) xH 


( 11 ) 


The above expressions in terms of the body-axis system assume 
a constant mass since the mass of the aircraft is not expected 
to change significantly in the short time during which data 
are collected. The components of the vector equations (10) 
and (11) are written as: 

F x = m (u +qw-rv) 

F y = m (v+ru -pw) (12) 

F z =m (w+pv- qu) 
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L=pI x -qI xy -±I xz +qr (I z ~J y ) + (r 2 -g 2 ) I yz -pql^rpl xy 
» = -pl^ql y -rl y ^rp (J x -J z ) * (p 2 -r 2 ) I^-qrt^+pqi: yz (13) 
W = -P J x z-^yz+^z+PZ <V J X ) + (<7 2 -P 2 > J *y- r P J yz +< ? rJ « 


B. EULER ANGLES 

The development in the previous section led to a set of 
six equations of motion derived for a reference system that is 
fixed to the airframe. Since the position and orientation of 
the airframe cannot be described relative to its moving 
reference frame, it is necessary to develop a way of defining 
the aircraft attitude relative to the earth. The angles that 
define the rotations that transform the orientation of the 
body frame to the earth-fixed inertial frame are called the 
Euler angles [Ref. 22:p. 88]. 

The transformation begins by aligning the body-axis system 
with the inertial system (Figure 2) . Looking down the Z axis 
into the X 1 Y 1 plane, a vector R in that plane can be resolved 
into x, and y, components. 

R=x^St + y,p (14) 

Rotate an angle Y about the Z axis to define a new coordinate 
plane X 2 Y 2 . The same resultant vector R can be expressed in 
terms of the components of the new coordinate system. 


20 







J?=X 2 * «■ y 2 f 


( 15 ) 


where 


x 2 =x 1 cosY+y 1 sinY 

y 2 =y 1 cos7-x 1 sinT (16) 


In matrix form the new coordinate system can be expressed as 
a product of the transformation matrix T* and the vector 
containing components x 1# y, and z 1# 


X 2 


cosT sinT 0 


*i' 

y 2 

= 

-sinY cosY 0 


Yy 

*2 


0 ° 1 


z y 


With the coordinate system now aligned with the X 2 , Y 2 , Z 2 
coordinate system, the same procedure is applied by looking 
down the Y 2 axis into the X 2 Z 2 plane and observing the R vector 
in terms of components x 2 and z 2 . Rotate an angle 0 about the 
Y 2 axis to form a new coordinate system X 3 , Y 3 , Z 3 . The vector 
R can now be expressed in terms of components x 3 , y 3 , z 3 by: 

R=x 3 ± + z 3 2 ( 18 ) 
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Figure 2 Euler Angle Rotation About the Z axis 

where in matrix form the new vector is a product of the 
transformation matrix T e and the vector in the x 2 , y 2 , z 2 
system. 

x 3 [cos© 0 -sin8] x 2 

y 3 - o l o y 2 (19) 

z 3 [sin8 0 cos8 z 2 

A similar procedure is applied to the X 3 , Y 3 , Z 3 coordinate 
system by rotating an angle * about the X 3 axis. The 







transformation matrix T # is found to be: 


10 0 


T* = 


0 cos® sin® 
0 -sin® cos® 


( 20 ) 


A complete transformation from the body axes to the 
inertial axes by three successive rotations about the Euler 
angles 7, 8 and • leads to the expression: 

X b =^T 6 T w X 1 ( 21 ) 

where X b is a vector of components in the body-axes system and 
X,- is a vector of components in the inertial-axes system. 
Figure 3 shows the earth-fixed axis relative to the body- axis 
system with the three successive rotations about the angles 7, 
e and t. 

Using the relationship in equation 3.21, the inertial 
velocities can be expressed in terms of the Euler angles and 
the components of velocity in the body-axis frame. The 
shorthand notation of References 21 and 22 where C e represents 
Cos 0 gives: 



( 22 ) 
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Figure 3 Euler Angle Orientation 

The angular rotation rate vector fl has body-axis 
components p, q and r where p is the roll rate about 
the longitudinal axis, q is the pitch rate about the lateral 
axis and r is the yaw rate about the vertical axis. Using the 
relationships of the transformation matrices T # , T e and T t 
given above, the body-axis components of the angular rotation 
rate vector A can be expressed in terms of the Euler-angle 
rates by: 






JP 


1 o -s e 


* 

q 

= 

0 Cq s*c 9 



r 


0 -S 9 C&Cg 


¥ 


( 23 ) 


Inverting equation (23) gives the Euler-angle rates in terms 
of the body-axis rotation rates: 


* 


i SinQTanQ CosQTarfi 


P 

d 

= 

0 CosQ -Sin<t> 


q 

¥ 


0 SinQSecQ CosQSecQ 


r 


(24) 


C. FORCE AND MOMENT EQUATIONS 

The applied forces on the aircraft can be expressed in 
terms of aerodynamic, gravitational and thrust components. 
The forces acting on the body of the aircraft relative to the 
earth are: 

F x = qsC x ~ mgSinB + T 

F y = qsC y +mgSin<l>Co& (25) 

F z = qsC z + mgCosQ CosB 

where q-bar is the dynamic pressure, s is the reference area, 
C x , C y , C z are the force coefficients relative to the aircraft 
state and T is thrust. 

The gyroscopic moment due to rotating parts in the engine 
is ignored and the aerodynamic moments are given by: 
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M x = gscC 2 

M y = qscC m (26) 

M 2 = qscC n 

where C t , C m and C n are non-dimensional coefficients of roll, 
pitch and yaw moments. 

D. STABILITY AND WIND AXES 

For stability and control analysis it is more convenient 
to consider the velocity of the aircraft relative to the air 
[Ref 18:p 8], Also it is more practical to express equations 
in terms of angle of attack, a and sideslip angle, B since 
these angles are measured directly during wind tunnel testing 
and flight testing. Angle of attack and sideslip are defined 
by: 

a * tan' 1 p=sin* 1 -£ (27) 

where V is the total wind-relative velocity with magnitude 
given by: 

V=Ju 2 + v 2 + w 2 (28) 


From Figure 4 the component velocities can be expressed in 
terms of the angle rotations through alpha and beta and the 
total wind-relative velocity V by the following process. 
Rotate about the w axis by an angle B: 
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cp -sp Ol u i 
sp Cp 0 V, 


Rotate down about the v axis by an angle a: 


\Ca 0 -Sal 


[Sa 0 Ca I w- 


The component velocities can then be expressed by: 


v = T.r„ o 



Figure 4 Flow Angles Relative to Alpha and Beta 
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Through a similar application of the angle relationships 
the axial forces can be expressed in terms of the non- 
dimensional coefficients C 0 and C L and related to the stability 
axes and the wind axes by a transformation matrix. 


C t 


Set 

-Ca 

c. 

Cp 


-Ctt 

-Sa 

Cm. 


c 

^wind 

C y*tn*. 


cp -sp 

Cp 

SP cp 

Cy 


(33) 


E. BODY AXI8 EQUATIONS OF MOTION 

Constructing the body-axis, six-degree-of-freedom 
equations of motion begins by differentiating equations (27) 
and (28) to obtain expressions for the time derivatives of V, 
alpha and beta. Equating equations (12) and (25) yields: 

u= -^C -gw+rv-gsin8 + — 
m * m 

v-=— C - xu + pw+gsin<bcos& (34) 

w r 

w= -2® c z -pv + qu + gcos<I>cos8 


Substitute equation (31) for u, v, and w and equation (34) for 
the time derivatives of u, v, and w. Use equations (32) and 
(33) to substitute in for the force coefficients Fx, Fy and 
Fz. The construction is completed by substituting equation 
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(33) to substitute in for the force coefficients Fx, Fy and 
Fz. The construction is completed by substituting equation 
(26) into equation (13). The resulting nine equations can 
then be expressed as: 


+ g (cos$cos8sinacos(J ♦ sin$cos9sin(J 

T 

-sinQcosacosP) cosacosp 

m 


(35) 


d = 


gs 


C L + q- (tanP) (pcosa +rsina) 


inVcosp 

+ —-2— (cos$cos8cosa +sin8sina) - rsina 
VtosP mVcosp 


(36) 


P =JZ£c v +psina-rcosa +-?cosPsin®cos8 
K mv V K 

h jsinji (gcosasin8-gsinacos^cosO + — cosa) 
V w 


(37) 



$I xy -iI xz =qsbC 1 +qr (I y -J r ) + (g 2 -r 2 ) I yz 

(38) 


+ PqI X z~ r Plxy 

~pl xy +ql y -il yz =< jscc m +xp{l z -l x ) + (r 2 -p 2 ) J xz 

(39) 


+ q r i xy ~pqi y z 


r qI yz +rI z =qsbC n +pq{I x -I y ) + (p 2 -g 2 ) 

(40) 


+rpl yz -qrl xz 


4> =p + gtan8sin0 +rtan8cos® 

(41) 


8 = gcos4> -rsin$> 

(42) 
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’P = rcos$sec8 + gsinG>sec9 


(43) 


The nine equations listed above are valid for a rigid 
vehicle in a constant wind with the following assumptions. 

1) Flat, non-rotating earth 

2) Constant mass 

3) Thrust vector along the X axis 

4) Body axis referenced to the aircraft center of gravity 


Singularities occur at e = ± 90 degrees, 6 = ± 90 degrees and 
at zero velocity. 


F. AERODYNAMIC MODEL 

Equations (12) and (25) can be combined to give a 
representation of the forces acting on the aircraft. Changes 
in the forces from a given reference condition can be 
represented by a Taylor series expansion as a function of the 
motion variables. In the development of the aerodynamic 
model, it is assumed that any change in the longitudinal force 
is a function of angle of attack, pitch rate and longitudinal 
control. Changes in the longitudinal forces and moment can 
then be represented by: 


ou aw do „ 


A6, 


A Au + ^Aw+^-A# 
ou aw aw 

AM= A u + 4^ Aw+ ~ A 
ou aw aw 


dM » BM a i 

r A *< 


(44) 
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Likewise it is assumed that the lateral-directional 
aerodynamic forces and moments are functions only of sideslip 
angle, roll rate, yaw rate and lateral-directional control- 
surface changes. Expansions for AY, AL and AN can then be 
expressed in terms of the motion variables Av, Ap, Ap, Ar and 
A 5 by: 


dv dp dr do z z 

at &L . ^ dL A „ ^ dL A , . dL A . . dL A * 

At at 4v * Tp Ap+ Ti Ajr * -W t A8 '* IT . A8 ‘ 

AW=i^Av+-^Ap + -^Ax + -^ AS + -^A6 
dv dp P dr db r r d6 a ‘ 


(45) 


In equations (44) and (45) a conventional three-axis 
control system consisting of aileron, rudder and elevator is 
assumed. Additional control terms may be added for non- 
conventional configurations. Reference 23, Chapter IV 
clarifies the significance of each of the derivative terms in 
the expansions and explains the physical significance by 
relating the longitudinal and lateral-directional forces and 
moments to the aerodynamic derivatives and motion variables. 

Stability and control derivatives in dimensional form are 
convenient for application of modern control theory methods 
when solving flight dynamic problems. References 21 and 22 
use the dimensional form in the development of the state-space 
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model. The non-dimensional coefficients, however, provide 
a convenient means of accounting for aircraft speed and size, 
and for air density since they are insensitive to dynamic 
pressure over large ranges [Ref. 24:p. 113]. Table 1 
summarizes the longitudinal stability and control derivatives 
in terms of non-dimensional coefficients and Table 2 
summarizes the lateral-directional stability and control 
derivatives. The derivatives represent the change in an 
aerodynamic force or moment coefficient due to a change in 
state or control input. 
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TABLE 1 LONGITUDINAL STATES AND DERIVATIVES 


Longitudinal Stability and Control Derivatives 

State or Control 

Normal Force 

Pitch Moment 

Angle of Attack 

CL a 

Cm ff 

Pitch Rate 

CLq 

Cm q 

Elevator Control 

CL 4e 



TABLE 2 LATERAL-DIRECTIONAL STATES AND DERIVATIVES 


Lateral-Directional Derivatives 


State or Control 

Side Force 

Roll Moment 

Yaw Moment 

Sideslip 

cy b 

ci B 

Cn B 

Roll Rate 

CYp 

C1 P 

Cn P 

Yaw Rate 

CY r 

ci r 

Cn r 

Aileron Control 


ci f . 

Cn «. 

Rudder Control 

C * 4r 

Cl,r 

Cn ir 
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IV. MAXIMUM LIKELIHOOD ESTIMATION 


A. ESTIMATION THEORY 

The problem of obtaining unknown aircraft stability and 
control derivatives through parameter estimation techniques is 
theoretically straight forward. In the previous chapter it 
was emphasized that the most important part of this analysis 
is deriving a reasonable mathematical model that consists of 
a set of dynamic equations that contains the unknown 
parameters. The system can then be excited by a pre¬ 
determined input specifically designed so that the response of 
the system will allow the values of the parameters to be 
determined. The system response to the pre-determined input 
is measured. Values of the parameters are then determined by 
the constraint that the computed response of the model match 
the measured response of the system. This deterministic 
approach to parameter identification is complicated by the 
presence of other elements such as process noise and 
measurement noise that must be considered in real-world 
systems. 

The mathematical model derived in the previous chapter is 
based on the physical laws that govern the system and is 
considered to be well defined. In deriving any model of a 
physical system it is necessary to make some compromise 
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between simplicity and accuracy. A model that is too 

complicated will require a higher expenditure of resources and 

a model that is too simple will not have sufficient accuracy 

to be useful. The objective is to define a mathematical model 

that represents the essential aspects of the physical system 

and will provide reasonably accurate predictions of the 

dynamic behavior [Ref. 21 :p. 67]. The term ’'reasonably 

accurate" implies that modeling errors are inevitable. In 

Reference 17 Maine and Iliff state: 

There is no comprehensive theory of modeling error 
available. Any modeling error is simply treated as state 
and/or measurement noise, in spite of the fact that the 
modeling error may be deterministic rather than random. 
This procedure has not been rigorously justified, but 
combined with careful choice of the model, is probably the 
best approach available. 

These unavoidable modeling errors are intentionally accepted 
in an effort to produce a simplified but practically useful 
simulation. 

The presence of both measurement noise and state noise 
changes the nature of the problem from a deterministic process 
to a stochastic process. Measurement noise arises from the 
realization that the system response cannot be perfectly 
measured and sensor inaccuracies are certain. State noise is 
random excitation of the system from external, unmeasured 
sources. The noise is a disturbance that acts uninterruptedly 
upon the system making it impossible to predict exactly how 
the system will behave. Consequently the unknown parameters 
must be estimated by some statistical criterion [Ref 17:p 8]. 
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The real value of the MMLE3 program is that the algorithm 
accounts for both state and measurement noise. 


B. MODIFIED MAXIMUM LIKELIHOOD VERSION THREE 
1. General System Model 

In the Modified Maximum Likelihood Estimation (MMLE3) 
method it is assumed that for the correct values of the 
unknown parameters, the system is correctly described by the 
dynamic model developed in the previous chapter. The model 
can be expressed in state-space format as a function of the 
states, inputs and unknown parameters represented by x, u and 
i respectively. The discretized model in state-space form is: 

= A k *it+B k u k +<* k ( 46 ) 

where the plant matrix A and the control matrix B contain 
values of the unknown parameters and u is the state noise 
vector. The initial condition is given by: 


Measurements are assumed to be made at discrete time 
intervals. The measurement equation is: 


z k = Cx k^k 


(48) 
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The matrix 

C is 

the 

measurement matrix and v 

k is 

the 

measurement 

noise 

vector. 

The state 

noise, e, 

and 

the 

measurement 

noise, 

v k* 

are 

assumed to 

be zero 

mean 

and 


uncorrelated [Ref. 25:p. 3]. In the presence of white noise 
the estimation is a stochastic process. For each estimate of 
the unknown parameters in the vector £, there is associated 
with the estimate a probability that the values in the 
response history are close to the measured values. The 
maximum likelihood estimates are the estimates that maximize 
that probability. 

2. Kalman Estimator 

A Kalman estimator is used to provide estimates for 
the states while minimizing the effects of the process and 
measurement noise. The difference between the measured value 
and the Kalman estimate at a given time, t(i), is called the 
residual and is given by: 




(49) 


The subscript £ indicates that the estimate is a function of 
the unknown parameters. To maximize the probability that the 
computed time history is close to the measured time history it 
is necessary to minimize the criterion given by [Ref 25:p 3]: 


J (l) -V c i>] r R ' 1 -V c i>] + f ln|*l (5°) 
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Where R is the innovations co-variance matrix. From equation 
(50) it can be seen that the effect of the state and 
measurement noise is only to increase the value associated 
with the criterion J(£). 

The discrete, stochastic estimation is a recursive 
process that makes an estimate after each measurement. The 
new estimate depends upon the prior estimate and the new 
measurement. In terms of a Kalman estimator equation, the 
prediction and correction steps can be written as: 


Xfc — + By, IZjj 

*k = 2 k + k[z k -z k ] 


(51) 


Equation (51) states simply that the corrected value is equal 
to the old estimate plus a confidence factor times the 
residual. The gain or confidence factor is chosen by solution 
of the Riccati equation such that the errors associated with 
the states are driven towards the convergence criterion as 
quickly as possible. MMLE3 uses only the steady state 
solution of the Riccati equation. A complete explanation of 
the discrete, stochastic process can be found in Chapter 6 of 
Reference 27. 

3. Cramer-Rao Bounds 

In any stochastic process some measure of accuracy 
is necessary. In the MMLE3 estimator the Cramer-Rao bound is 
a measure of the reliability of the estimate based on the 
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information obtained from the dynamic maneuver. In simpler 
terms, the Cramer-Rao bound is an estimate of the standard 
deviation of the estimates. The measure of reliability is 
improved when a sufficient number of measurements are provided 
and the assumptions used in the prediction model are valid. 
The Cramer-Rao bounds are the square roots of the diagonal 
elements of the information matrix. A precise mathematical 
description of the information matrix and the calculation of 
the Cramer-Rao bound can be found in References 25 and 27. 
For the purposes of this study it is enough to know that the 
Cramer-Rao bound is a good indicator of the accuracy of an 
estimated parameter. 

4. A-Priori Information 

Wind tunnel data from the study conducted by the 
NASA Langley Research Center, Flight Dynamics Branch, provide 
some insight into the aerodynamic characteristics of the BQM- 
147A. Measurements tabulated in Reference 2 are useful in 
providing estimates for some of the stability and control 
parameters. The MMLE3 parameter estimation routine allows the 
estimation algorithm to consider this a-priori information by 
adding a quadratic penalty function to the criterion [Ref. 
27:p. SR-4]. Equation (50) can then be re-expressed by: 
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( 52 ) 


N 

JU) «—-2 € (t.)] T R" 1 [*(t f ) -f<(t.)] + 

2 i«i 

i(5-c 0 ) T w a-i 0 ) 

2 

where £ 0 is the vector of a-priori values and W is the 
weighting matrix. The a-priori information can be used to 
speed the convergence; however, the results will be biased 
towards the weighted estimates [Ref 27:p ST-8], A good 
technique to avoid the problems associated with biased 
estimates is to remove the weighting after initial 
convergence. The initial weighted values should be close 
enough to the maximum likelihood estimates to allow rapid 
convergence to an unbiased estimate when the weighting is 
removed and additional iterations are allowed. 

5. Linearized Equations of Motion 

In the analysis of aircraft stability and control 
it is often more practical to restrict the area of interest to 
steady-state motion and the response of the aircraft to small 
perturbations about a steady-state condition [Ref 24:p 33]. 
Equations (35) through (43) in Chapter 3 are a full set of 
coupled, non-linear equations that completely describe the 
dynamic behavior of the aircraft. The simplification of the 
model by applying small-angle theory significantly reduces the 
complexity of the model without necessarily compromising 
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accuracy. For a wide range of practical applications on well- 
behaved systems, a linear model is sufficient. 

References 22 through 25 provide details in applying 
small-angle theory to the non-linear equations of motion from 
which the linearized equations used in MMLE3 are derived. 
From equations (36), (39) and (42), a small angle 
approximation leads to the simplified longitudinal equations: 


4*- 



■q+& 
H V 


(53) 


4-&c m 


(54) 




(55) 


In the approximation, a plane of symmetry about the X-Z plane 
is assumed, making I yz and I equal to zero. The aerodynamic 
coefficients in the above equations are: 


~ + C^b + q. 


C =C m a+ C- b + CL +C_ 
m m m m g 2 V m< ‘ 


(56) 


(57) 


From equations (56) and (57) the vector £ is expressed as a 
function of the unknown stability and control derivatives by: 






5 = 5 (, C Lt , C,^, C,^, C a< ) 


( 58 ) 


From the application of the small angle approximation to 
equations (37), (38), (40) and (41), the simplified lateral- 
directional equations are: 

¥~ r < 59 > 

pI K -rl xz = qsbC 1 + qz (I y -I z ) +pql xz ( 60 ) 

tl z -pl xz = qsbC n + pq ( I x - I y ) - qz I xz ( 61 ) 


4>=p + rcos<|>tan0 + gsin<J>tan0 


(62) 


The aerodynamic coefficients for the lateral-directional 
equations are: 


c„=c„ Q + c, + c„ - r ^ > - + c, & + c, 

y r p2V r * 2V *» 1 


(63) 


C, = C, 6 + C, ^ + C, — + C, 6 + C, 
J J » K 2 i>2V lr 2V j » j 


(64) 
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(65) 


rb 


■C n fi + C„ & + C„ — + C_ b *C n 

n^r n p 2 v n r 2 V ^ 110 


The "delta" parameter is a general expression for control 
input and must be expanded to consider all inputs that 
contribute to the total force or moment coefficient. The 
contributions to the total side force coefficient from roll 
rate and yaw rate are small compared to the other contributing 
parameters and are usually ignored. The vector of unknown 
parameters, can be expressed as a function of the lateral- 
directional stability and control derivatives by: 


(C Y . C, 




C K i Ci i 
A P x x 


i 




( 66 ) 


C. pEst 

For many aircraft applications the linearized 
equations of motion provide a reasonably accurate model of the 
system for small disturbances. For some applications however, 
the linear model is not sufficient and a more accurate non¬ 
linear model is required. Flight testing at the Dryden Flight 
Research Facility highlighted the limitations of the linear 
model when applied to maneuvers at extreme flight conditions. 
In 1987 Murray and Maine [Ref 9] introduced the pEst program 
for parameter estimation in non-linear dynamic systems. 
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l. General System Model 


Conceptually, parameter estimation in non-linear 
dynamic systems is the same as the parameter estimation in 
linear systems. The general system model is separated into a 
continuous-time state equation and a discrete response 
equation represented by [Ref. 9:p. 3]: 

±<t) =f [X(t) , U(t) , 5J (67) 

z(t d ) = g[x(t d ), u(t d ) , (68) 

The cost function is the criterion by which the model's 
computed response is measured against the measured time 
history response. The pEst program defines quantitatively the 
criterion by: 

i 

J(l) = — T [ Z ( fcj) -2(t d )) T W [z (t d ) -zU d )] (69) 

2 n z n t 


Where n z and n t are the number of response variables and the 
number of measured time history points. 

2. Estimator 

The feedback process, like that used in the MMLE3 
estimator, is recursive. The non-linear equations used in 
pEst preclude using the discrete-time Kalman estimator 
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formulation of MMLE3; however, the feedback term is still 
proportional to the difference between the measured response 
and the computed response [Ref. 9:p. 3]: 

St( t s ) =x( tj *k [z( tj) -2 (tj) ] (70) 

3. Parameters 

The parameters are the values that can be estimated by 
pEst. Each of the parameters has five attributes: current 
value, estimation status, predicted value, Cramer-Rao bound 
and change in value from the previous iteration [Ref 9:p 11]. 
The initial current value is a best estimate of the parameter 
entered by the user. Typically the initial current value 
input is derived from a source such as DATCOM, [Ref. 30 ] or 
it may be the result of measured data from either wind tunnel 
studies or flight tests. The estimation status indicates to 
the user which parameters are being estimated. A parameter is 
either active or inactive. The predicted value is the value 
of the parameter computed by pEst and it is used to update the 
current value for each successive iteration. The Cramer-Rao 
bound and the change in parameter value are defined in the 
estimation process. The Cramer-Rao bound is a measure of the 
accuracy of the estimate. A list of the names and a 
description of the parameters available in pEst can be found 
in Section B.16, Reference 9. For this study, the parameter 
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vector, i, is the same as the parameter vector used in the 
MMLE3 routine. 

4. States and Responses 

The attributes of the states are status and 
integration limit. The attributes of the response variables 
are status and weighting. If a state is active, the state 
equation is integrated and its value is used in the equations 
of motion. If the value of the state exceeds the integration 
limit, the integration is terminated. If the response 
variable, or output variable, is active then its time history 
values are computed and placed in a special output file. The 
weighting of a response variable specifies the value to be 
used in the W matrix of the cost function. Weighting the 
response variables based on a-priori knowledge may help the 
routine converge more rapidly; however, the results may be 
biased towards the weighted estimates. Table 4.13 of 
Reference 24 shows the relative importance and prediction 
accuracies of the stability derivatives using theoretical 
methods. The table may be used as a guideline if weighting is 
desired. 

5. Equations of Motion 

Equations (35) through (43) define the non-linear 
equations of motion used in the pEst routine. The assumption 
of symmetry about the X-Z plane that was used in MMLE3 is not 
used in pEst; however, the assumptions of fixed aircraft 
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geometry and constant mass are valid [Ref. 9:p. 17]. In this 
study, equations (36) through (42) are integrated to obtain 
the predicted state variables in the vector x-tilde. A list 
of the 14 response equations can be found in Section 6.1.4 of 
Reference 9. 
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V. MMLE3 APPLICATION 

Recent thesis work by LCDR Robert Graham [Ref. 6] 
incorporated a PC-based parameter estimation capability into 
the NPS research program. The Maine and IIiff MMLE3 code has 
been used for a number of years for large scale work where the 
use of mainframe computer systems was justified [Ref. 9:p. 1]. 
Research by Erickson [Ref. 8] showed that for less demanding 
work, a PC-based parameter estimation routine for a 
longitudinal model generally provides comparable results when 
compared to more complex mainframe programs. The code written 
by Graham uses the commercially available software package, 
MATLAB*, to build interactive files for both the longitudinal 
and lateral-directional models. The MATLAB State-Space 
Identification Tool, (SSID) [Ref. 28], comprises a main 
Modified Maximum Likelihood Estimator (MMLE) program, along 
with several supporting functions for parameter estimation. 
The supporting functions are used to convert the parameter 
vector into the state-space model and to arrange the necessary 
inputs for the MMLE program. A complete listing of the code 
written by Graham can be found in Appendices E and F of 
Reference 6. The modifications to the supporting functions in 
Reference 6 required for this application are discussed in the 
following sections. Figure 5 [Ref. 9:p. 43], shows how the 
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various MATLAB files interact in the parameter estimation 
routine. 



Figure 5 Interactive MATLAB .M Files 


A. DATA ACQUISITION 
1. Flight Test 

The data made available for this study were provided 
by NASA Langley Research Center and were the result of flight 
tests conducted at the NASA Langley Plumtree Test Site in 
March, 1992. This flight test comprised a single set of 
flight maneuvers consisting of pitch, roll and yaw doublets. 

A partial sequence of events is listed in Table 3. The roll 
and yaw doublets were merged into a single input file so that 
all of the lateral-directional modes would be excited. 
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Figures 6 and 7 show a typical control input for the doublet 
maneuvers. 

2. Atmospheric Data 

Specific atmospheric data for the date and time of the 
flight tests were unknown. Best estimates of atmospheric 
conditions for the test site were from a local area weather 
briefing and were recorded on the flight card by the flight 
test pilot. 

3. Sensor Placement 

Twin alpha and beta sensors were symmetrically located 
about the longitudinal axis and 2.5 feet in front of the 
vehicle center of gravity (eg). An averaging scheme was used 
on the two sets of measured data and a single set of data for 
alpha and beta was provided. Accelerometers for X, Y and Z 
axes were located on the eg point. A rate sensing gyro was 
located on the waterline, aft of the eg. Details on sensor 
characteristics were not available. 

The aircraft equations of motion developed in Chapter III 
assume that the sensors are located at the eg. Since several 
of the observed signals are functions of the position of the 
sensors, it is necessary to correct the sensor measurement for 
the offset. For the longitudinal model the measured variables 
a,,, and A n are corrected to the eg by [Ref. 18:p. 26]: 


a =a + a P q 
eg m „ M 


(71) 


50 





(72) 



-Ar 


^g-fL V 

9 9 


For this application, is positive forward of the 
eg and was given a fixed value of +2.5 feet. Since the X-and- 
Z axis accelerometers were located at the eg, the value of X,,, 
and Z,,, is zero. For the lateral-directional model, the 
measured variables 13^ and Ay are corrected to the eg by [Ref 
18:p 27): 



(73) 




(74) 


The X distance to the beta probe was set at +2.5 feet and the 
value of X_„ and Z av was set to zero. 

4. Measured signals 

A total of twenty-four signals were measured and 
recorded during the flight test. These signals are presented 
in Table 4. Details of how the data were filtered, digitized, 
time tagged and recorded were not available. 
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B. USER FUNCTION 

The PC-based MMLE3 application requires a user-written 
function to convert the parameter vector into a state-space 
system model [Ref. 28:p. SR-3]. The following two sub¬ 
sections describe the creation of the P2SS.m file for this 
study. 

l. Longitudinal Modal 

In the general system model given by equation (46), 
the response of the system is a function of the state vector, 
x, and the control input, u. Since the measured flight test 
data do not include pitch angle and in the short-period 
approximation forward velocity is assumed to be constant, the 
simplified longitudinal model can be reduce to equations (53) 
and (54) for a and q. These two equations contain terms which 
are neither functions of the state variables or the control 
inputs. Consequently, the control input was augmented with a 
unity input to incorporated these terms in the state-space 
model. The state vector is represented by: 



The control input is represented by: 



( 76 ) 
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In the state-space formulation the vector p is the vector 
of stability and control parameters. This vector is expressed 
in terms of its elements by: 

P = [ C L« ^m,,] (77) 


The plant matrix can now be expressed in terms of the elements 
of the parameter vector by: 



-rr^—-p(l) 

gwVcosctj 

lp p (2) 



(78) 


The control matrix is expressed in terms of the elements of 
the parameter vector by: 



-25—p (4) 

gwVcostty^ 


5) 

1 y 


QS 

gwV cosffj 
^ 3 . — (p (2 ) a x 

I y 


(p( 1) +p(4) 6 ej ) -g x 
+ P< 5 > 'fff~P (3> 


(79) 


In the longitudinal application it is desired to evaluate 
the dependence of any change in normal acceleration upon 
longitudinal control input. The equation for normal 
acceleration in terms of non-dimensional stability and control 
parameters is derived in Chapter 4 of Reference 22. The 
general system output equation can be expressed in terms of a 
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parameters is derived in Chapter 4 of Reference 22. The 
general system output equation can be expressed in terms of a 
parameter vector, y, that contains the variable A n in addition 
to the states a and q. The output and measurement equations 
of the form given in equation (48) are augmented by a feed¬ 
through matrix, D. The output matrix, C, can be expressed by: 


c= 


1 

0 


QS 

mg 


p( 1) 


o' 

1 

0 


(80) 


The feed-through matrix, D, is given by: 



0 

0 


-2^p(4) 

gw 



(81) 


The subscripted states, a, and q t , are the elements of the 
initial state vector x 0 . The input for the initial vector 
comes directly from the first element in the measured data 
file for each of the states. 

The longitudinal model given above specifies the state 
equation in continuous form. In computing the continuous-to- 
discrete-time conversion, the sample time, dt, is required. 
Since MATLAB will not allow a global variable to be used in a 
function argument, the P2SS.ro file must be modified to allow 
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2. Lateral-Directional Model 


In the lateral-directional model the state vector, x, 
is represented by: 


x= [6 p r 0] T 


(82) 


The control input vector, u, is represented by: 


u = 



(83) 


The subscripted variables S a and S r represent the aileron and 
rudder control inputs. The input vector is augmented by a 
unity input to incorporate the additional terms that influence 
lateral acceleration in the measurement equation. The 
measurement vector is: 

y=[fipr^A y ] T (84) 


The lateral-directional parameter vector, p, is: 


P-[ C r. c ,. e „. c i 


s c ‘. 


c 







(85) 


The only change to the NPSP2SS4.M file in Reference 6 


55 




was to input the time step, dt, for the sampling rate of the 
measured data file. 

From Table 4 it can be seen that roll-angle was not 
measured during flight test. Since the spiral mode is 
dominated by the roll angle component, it was necessary to 
numerically integrate the roll rate, p, in order to provide a 
roll angle input estimate. For N intervals with equal 
spacing, the extended trapezoidal rule can be written as: 




N-1 


f(a) +2£ f ( a+ J ( d t)) (*>) 


*E 


( 86 ) 


A complete development of the extended trapezoidal rule and a 
discussion of the error term, E, can be found in Chapter IV of 
Reference 29. 

C. INITIALIZATION 

The initialization macro is a supporting function that 
imports the measured time history file, prompts the user for 
the physical dimensions of the aircraft and prompts for the 
initial values of the parameters that are to be identified, 
l. Physical Characteristics 

Table 5 presents the aircraft physical characteristics 
that were input for this study. To be consistent with the 
physical characteristics used in the NASA static wind-tunnel 
tests, the theoretical root chord was used as the reference 
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1. Physical Characteristics 

Table 5 presents the aircraft physical characteristics 
that were input for this study. To be consistent with the 
physical characteristics used in the NASA static wind-tunnel 
tests, the theoretical root chord was used as the reference 
wing chord. Moments of inertia and sensor offset information 
were provided by the Marine Corps UAV Project Office, 
Intelligence Systems, Quantico, VA. 

2. Initial Estimates 

The SSID Modified Maximum Likelihood Estimator 
requires a set of predicted aerodynamic characteristics to be 
used as starting values for the iterative algorithms. In this 
study a combination of wind tunnel data and analytic estimates 
was used as a source of predicted values. Tables 6 and 7 
present the parameter, the predicted value of the parameter 
and the source of the predicted value for the longitudinal and 
lateral-directional case studies. 

D. MODIFIED MAXIMUM LIKELIHOOD ESTIMATOR 

The Modified Maximum Likelihood Estimator is the main 

program used by the MATLAB SSID tool and it is the only 

program that is called by the user. MMLE calls the supporting 

functions to provide the necessary inputs for the estimation 

algorithm. The primary inputs used by the MMLE macro are: 

1) uydata: An n-by-m matrix where the m columns are the 
measured time history inputs of the vector u followed by the 
outputs of the vector y and n is the number of time steps. 
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2) p2snam: The name of the P2SS function that converts the 
parameter vector to the state-space model. 

3) pO: The initial parameter estimates entered during 

initialization. 

4) rmsO: A vector containing the standard deviations of a- 
priori parameter estimates. 

5) ggO: The initial guess of the innovations co-variance 
matrix. 

6) linesearch: An option to help convergence during the 
Constrained-Newton minimization phase. 

7) opt: A vector that specifies the number of iterations 
for the various minimization phases and the convergence 
criteria that permits termination of each phase. 

For this study, rmsO was a zero vector indicating that 
initial parameter weighting was not used. The initial 
innovations co-variance matrix, ggO, was set to a small, non¬ 
zero diagonal matrix to minimize the effects of the noise. 

There are four variations of the Newton method available 
for the minimization of the weighted error sum. The opt 
vector was set so that minimization begins with a Marquardt, 
gradient algorithm. The gradient algorithm performs well when 
the cost function is very steep, but performance tends to 
deteriorate as a minimum is approached and the cost becomes 
more flat [Ref. 28:p. SR-8]. When the user-specified narq 
value is reached, the process has been conditioned well enough 
for a Constrained-Newton minimization algorithm to begin. 
This algorithm is most accurate near the minimum of a cost 
function when the cost function is quadratic [Ref. 9:p. 16]. 
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The linesearch option is set "true" to help convergence in the 
Constrained-Newton phase. Reference 28 discusses each of the 
minimization phases in more detail. For this study the value 
of marq was set at 0.02 and the convergence criteria was set 
at 0.001. All parameters are identified in both the Marquardt 
and the Constrained-Newton phases. 
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TABLE 3 FLIGHT 03-24 FLIGHT CARD 


Pre-Flight 

and Takeoff 

P01 

Engine Off Data 

P02 

Static Measurement (nose up) 

P03 

Static Measurement (right wing up) 

P2A 

Takeoff Data 

Stability and Control Maneuvers (100%) 

P2B1 

Pitch-Up Doublet 

P2B2 

Pitch-Down Doublet 

P2B3 

Roll-Right Doublet 

P2B4 

Roll-Left Doublet 

P2B5 

Yaw-Right Doublet 

P2B6 

Yaw-Left Doublet 

Trim Flight 

Data 

P2E1 

20% Throttle 

P2E2 

50% Throttle 

P2E3 

100% Throttle 
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TABLE 4 FLIGHT 03-24 MEASURED SIGNALS 


Columns 1 through 4 


ALPHA 

BETA 

VELOCITY 

RPM 

Columns 5 through 8 

RUDDER 

AIL_LEFT 

ELEV_LEFT 

AILRIGHT 

Columns 9 through 12 



ELEV_RIGHT 

THROTTLE 

PITCH_RATE 

ROLLRATE 

Columns 13 through 16 



YAW_RATE 

AX 

AY 

A2 

Columns 17 through 20 



AIL_CMD 

RUDDER_CMD 

ELEV_CMD 

ALTITUDE 

Columns 21 through 24 



SERVO VOLT 


SERVO AMP 


BLANK 


BLANK 



























































TABLE 5 BQM-147 PHYSICAL CHARACTERISTICS 


Reference Area (S) 

21.24 ft 2 

Reference Chord (c-bar) 

4.458 ft 

Wingspan (b-bar) 

8.167 ft 

Gross Weight (w) 

65 lb 

Moment of Inertia (Iy) 

3.443 slug-ft 2 

Moment of Inertia (lx) 

3.648 slug-ft 2 

Moment of Inertia (Ixz) 

0 

Moment of Inertia (Iz) 

5.965 slug-ft 2 

Sensor Offset (Xap) 

+2.5 ft 

Sensor Offset (Xbp) 

+2.5 ft 

Sensor Offset (Zan) 

0 

Sensor Offset (Zay) 

0 

Sensor Offset (Xan) 

0 

Sensor Offset (Xay) 

0 
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TABLE 6 INITIAL PARAMETER ESTIMATES, LONGITUDINAL 


Parameter 

Value 

Source | 

C L. 

2.4350 

Wind Tunnel [Ref 2:p 20] 

Cm a 

1.6900 

Wind Tunnel [Ref 2:p 20] | 

Cm q 

-0.1089 

DATCOM [Ref 30] 1 

C L«e 

-0.6670 

DATCOM [Ref 30] I 


-0.4060 

DATCOM [Ref 30] 
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TABLE 7 INITIAL PARAMETER ESTIMATES, LATERAL-DIRECTIONAL 


Parameter 

Value 

Source 

c ro 

-.2005 /rad 

Wind Tunnel [Ref 2:p 24] 

C IB 

-.0325 /rad 

Smetana [Ref 31:p 108] 

C„ B 

.0430 /rad 

Wind Tunnel [Ref 2:p 24] 

C t P 

-.3525 /rad 

Smetana [Ref 31:p 122] 

C np 

.0027 /rad 

Smetana [Ref 31:p 122] 

c lr 

.0729 /rad 

Smetana [Ref 31:p 133] 

C nr 

-.0230 /rad 

Smetana [Ref 31:p 135] 

C Ua 

.1377 /rad 

Smetana [Ref 31:p 124] 

Cr^a 

-.0194 /rad 

Smetana [Ref 31:p 136] 

C Y*r 

.1130 /rad 

Smetana [Ref 31:p 145] 

C lSr 

.0023 /rad 

Smetana [Ref 31:p 147] 

Crfr 

-.0138 /rad 

Smetana [Ref 31:p 149] 
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VI. pEst APPLICATION 

For many aircraft applications the linearized equations of 
motion provide a reasonably accurate model of the system for 
small disturbances about an equilibrium condition. For some 
applications, however, the linear model is insufficient and a 
more accurate non-linear model is required. Flight testing at 
the Dryden Flight Research Facility highlighted the 
limitations of the linear model when applied to maneuvers at 
extreme flight conditions and with some unique aircraft 
configurations that exhibited non-linear dynamic behavior. In 
1987 Murray and Maine introduced a program for parameter 
estimation, pEst, in non-linear dynamic systems. A full 
description of the pEst program can be found in Reference 9. 

The pEst program interfaces directly with two other 
programs. The GetData program, Reference 32, is a FORTRAN 
utility program that is used to extract selected signals from 
a time-history input file and to write the selected signals in 
a user specified format to an output file. The second program 
interfacing with pEst is XPLOT, Reference 33. The XPLOT 
program is a plotting utility that is used to plot selected 
signals from the measured time-history file and the computed 
time-history file. The following sections describe how the 
pEst program was set up for the estimation of the BQM-147 
stability and control derivatives. 
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A. MEASURED TIME-HISTORY DATA 

The twenty-four signals that comprise the recorded time- 
history data from the flight test, Table 4, were re-formatted 
so that they could be read into the GetData program. The 
FORTRAN utility used for re-formatting is listed in Table 8. 


TABLE 8 DATA CONVERSION PROGRAM 


program data 
dimension dat(20) 

open(l,file='data01',status=•old 1 ) 
k=0. 

read(1,*) 

20 read(l,30,end=97) (dat (i),i=l,24) 

30 format(f6.0,2f7.0,f8.0,5f7.0,3x,f4.0,9f7.0, 

# Ix,f7.0,f6.0,3f7.0) 
t=k/16. 

write(2,60) t,(dat(j),j=l,24) 

60 format(f10.3,lOx,3g20.14/5(4g20.14/),g20.14) 

k=k+l 
goto 20 

97 write(6,120) 

120 format(///' **** Normal Termination ***'///) 

99 stop 

end 


The GetData program was used to select the desired signals 
from the available time-history data. The signals not 
applicable to the pEst program were deleted. Those signals 
specifically required for the computation of the estimates 
were renamed so that the signal names on file matched the 
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signal names expected by pEst. One new signal was added to 
provide an input for dynamic pressure. The GetData output 
file was saved in a compressed format for use by the pEst 
program. 

B. OPERATIONAL STATUS 

The pEst program uses a file named current to store the 
operational status of the program. The current file is used 
to provide initial values for program variables, to store the 
status of the program for recovery and to store the results at 
the end of each run. The primary elements of the current file 
are parameters, states, outputs and constants. A listing of 
the current parameter status for the longitudinal and lateral- 
directional applications can be found in Tables 9 and 10. 

1. Parameters 

The parameters are variables that can be estimated by 
pEst. The five attributes of each parameter are: 

1) Current value: Typically entered by the user 

2) Predicted Value: Computed by pEst 

3) Estimation Status: Active (T) or Inactive (F) 

4) Cramer-Rao bound 

5) Change in value from previous iteration 

The current values of active parameters were the same 
as those used in the MMLE3 routine multiplied by a constant to 
reflect the change in units required for use in pEst. The 
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active parameters were those that are listed as elements of 
the parameter vectors given by equations (77) and (85). Only 
those parameters with an estimation status of "true" are 
estimated by the program. 

The parameter list also allows for a sensor placement 
input. In the MMLE3 routine the sensor offset was corrected 
to give a signal that simulated a measurement at the eg. In 
pEst, the effects of the sensor position are included in the 
observation equations. A list of the observation equations 
can be found on page 20 of Reference 9. A description of the 
available parameters can be found in Section B of Reference 9. 

2. States 

The status of the state variables in the current file 
is used to determine which of the state equations are to be 
integrated. The state equations used by pEst are equations 
(35) through (43) . The state vectors for the longitudinal 
model and the lateral-directional model are given by equations 
(75) and (82) . Since roll angle was not measured, the state, 
phi, was set inactive for the lateral-directional application. 

3. Outputs 

The attributes of the outputs are status and 
weighting. Default weighting was used in the initial current 
file in order to avoid biasing the solution in a particular 
direction. The outputs are the elements of the output vector, 
y, of the form given by equation (48). 
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4. Constants 


The constants are variables that cannot be estimated 
by pEst. In this application, constants were used to input 
the physical characteristic> -•'f the BQM-147. The constants 
and their respective values are the same as those listed in 
Table 5 for the MMLE3 application. 

C. MINIMIZATION ALGORITHMS 

The pEst program has several minimization algorithms 
available. In this application a default Nevrton-Raphson 
minimization routine was used. After each iteration the 
active parameters, the current values and the change in 
current value from the previous iteration are displayed. The 
Cramer-Rao bound is computed with the Newton-Raphson algorithm 
only. The Cramer-Rao bounds are displayed at the user's 
discretion. 

D. COMPUTED TIME-HISTORY 

For each output whose status is active a variable time- 
history is computed. The 'write' command is used to write 
the time history into a user specified output file. The 
signals written to the output file are the computed states, 
the computed responses and the residuals. The output file can 
be read by the XPLOT program and user selected signals can 
then be called for display. 
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TABLE 9 pEst INITIAL STATES; LONGITUDINAL 


File Name: Current.p2bl j 

Active Parameters 

Value 

cNormO 

1.637171976 

cNorma 

0.046710000 

cNormde 

0.029500000 

cmO 

0.0 

cma 

-0.001900000 

cmq 

-0.677000000 

cmde 

-0.007100000 

alphaO 

0.0 

qO 

0.0 

ka 

1.0 

Constants 

Value 

mass 

2.0300 

ix 

3.6480 ' 

iy 

3.4430 

iz 

5.9650 

ixy 

0.0 

iyz 

0.0 

ixz 

0.0 

area 

21.2400 

span 

8.1670 

chord 

4.4580 

States 

Status 

V 

F 

alpha 

T 

q 

T 

theta 

F 

Outputs 

Status 

V 

F 

alpha 

T 

q 

T 

theta 

F 

an 

T 
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TABLE 10 pEst INITIAL STATUS; LAT-DIR 


File Name: Current.uav 


Active Parameter 

Status 

cyO 

0.0 

cyb 

-0.00349900 

cydr 

0.00197210 

clO 

0.0 

clb 

-0.00056719 

clp 

-0.35248000 

clr 

0.07290000 

clda 

0.00240300 

cldr 

0.00007888 

cnO 

0.0 

cnb 

0.00075044 

cnp 

-0.00270000 

cnr 

-0.02300000 

cnda 

0.00033857 

cndr 

-0.00055497 j 

kb 

1.0 

States 

Stati -f 

beta 

T 

P 

T 

r 

T 

phi 

T 

Outputs 

Status 

beta 

T 

P 

T 

r 

T 

phi 

T 

ay 

T 

Note: Constants are the 

same for both 

| Longitudinal and Lateral-Directional | 
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VII. RESULTS AMD DISCUSSION 

A stated purpose of this study was the practical 
application of parameter estimation for estimating stability 
and control derivatives from BQM-147 flight test data. 
Previous information from a wind-tunnel study was made 
available; however, the results of that study did not provide 
estimates of all the required parameters, and changes in the 
production configuration would be expected to have some 
influence on the value of the parameters. Consequently, 
parameter estimation was not used to validate wind-tunnel or 
analytical predictions because there was no way to look for 
consistency between the data sets, nor was there sufficient 
flight-test data for a valid comparison. The complete 
aerodynamic data-base was required to be built using flight 
test data only. The wind-tunnel and analytical predictions 
were used only to provide starting values for the iterative 
algorithms and to provide a reasonable range for flight test 
results. 

When a complete data base must be built using flight data 
only, extensive coverage of the flight envelope is needed. 
The flight test data made available for this study comprised 
a single set of flight maneuvers consisting of pitch, roll and 
yaw doublets. For the longitudinal case, the short-period 
mode was excited by an elevator doublet. Table 3 shows that 
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data from only two maneuvers, a pitch-up doublet (P2B1) and a 
pitch-down doublet (P2B2), were available for exciting the 
short-period mode. The lateral-directional model included the 
roll, dutch-roll and spiral modes. The best way to excite all 
three of these modes is to have both aileron and rudder 
control inputs. Since the lateral-directional maneuvers were 
performed under similar flight conditions, the individual roll 
(P2B3) and yaw (P2B4) doublets were paired into a single 
maneuver. The measured data describing the maneuvers and the 
responses are found in Appendix B. 

A. MMLE3 

The PC-based MMLE3 routine used in this application relies 
on the locally linearized aerodynamic model presented in 
Chapter V. The results of the wind-tunnel study [Ref. 2] show 
that the behavior of the vehicle is approximately linear for 
angles of attack up to about 15 degrees. The linearized model 
is valid only when the assumptions made in the linearization 
process are valid. Thus, the model is ideally suited for 
small angle perturbation maneuvers and it can be expected that 
as the magnitude of the maneuver increases, the locally 
linearized model loses its validity and the estimates from the 
model are a less accurate representation of the aircraft 
dynamics. 

The process noise option was used in the MMLE3 routine. 
The risk in using the process noise option is that the Kalman 
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gains tend to become high, giving the high-frequency modeling 
errors more opportunity to influence parameter values. 

1. Model Validation 

In Chapter III, the changes that were necessary to 
Graham's routine [Ref. 6] to reduce the order of the model 
were discussed. In order to validate the changes, the 
modified model was run using the same F-14A data that were 
presented in Reference 6. The results of the run shown in 
Figures 8 through 11 compare well with Graham's results. 

2. Longitudinal Model 

The computations in the MMLE3 routine were limited by 
the working precision of MATLAB. Prior experience showed that 
the PC-based routine was highly sensitive to the initial input 
derivatives. After several adjustments to the initial 
parameter vector, the run using the P2B1 maneuver converged 
successfully. The run using the P2B2 maneuver did not 
converge because of a singularity due to working precision. 
The results of the P2B1 maneuver are shown in Figures 12 
through 15. The MATLAB output listing showing the input and 
determined derivatives, per radian, is found in Table 11. 

The estimated response for a appears to lag the measured 
response by a time step at the beginning of the maneuver and 
is unable to keep up with the amplitude of the measured 
response in the pitch-down portion of the maneuver. This 
behavior could have been predicted since pitch angle, which 
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was not measured, is an important term in the a state equation 
when the maneuver is not restricted to small perturbations 
from a stable flight condition. The estimated q response 
generally correlates well with the measured q response though 
it does not pick up the high rates at the extreme ends of the 
maneuver. The differences in q may be due to the ignored 
terms in the state equation in the linearization process. The 
estimated normal acceleration response appears to correlate 
well. 

3. Lateral-Directional Model 

Inspection of the data in Appendix B indicates that 
the data were exceptionally noisy in the absence of any 
control input. Noise in the general frequency range of the 
system response was a problem. When the P2B3 and P2B5 
maneuvers were paired for the lateral-directional run, the 
process noise was too large and the Gauss-Newton iteration was 
terminated. Consequently there were no results for the 
lateral-directional model in the MMLE3 routine. 
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Figure 12 BQM-147 Elevator Input (P2B1) 



BOM 147 ESTIMATED AND MEASURED AOA RESPONSE 

15 

_ - -r - , -- » • • • ‘ - - * - * - r -, - • - - 



• * Measured Data Points 


in 

• " „ o Estimated Response 



ft On 


5 

ynnnnoyoy900*onnOOn° f ° # n n o n nn n 0 n nr* n n non n g o Q O O O * f ' 

• * * * . . « • 

# •• *••**••••• 


I 

n • % 

n 


Hi 

3 n 

• 

n 


< 



Q 

n 


< 

n 


-5 

• • 


-in 

• 



._t_ J _ . _ x .. - x - • - — ■ * - 


n 0.5 1 1.5 2 2.5 3 3.5 4 4.5 j 


Time (seconds) 


Figure 13 BQM-147 a Response (P2B1) 
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Figure 15 BQM-147 Normal Acceleration Response (P2B1) 
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TABLE 11 MATLAB RESULTS FOR P2B1 


NPS PARAMETER IDENTIFICATION MACRO FILE 
FOR ACTUAL FLIGHT TESTS USING SIMPLFIEC SHORT PERIOD 
LONGITUDINAL STABILITY AND CONTROL DERIVATIVES 








RESULTS 

pid 

p(pid) 

pref 

cramer 

2fcramer 

insens 

gdop 

1.0 

2.7391 

2.4350 

0.1084 

0.3085 

0.0780 

1.3902 

2.0 

-0.0887 

-0.1089 

0.0028 

0.0079 

0.0018 

1.5088 

3.0 

-0.1078 

-0.6770 

0.0414 

0.1178 

0.0159 

2.6088 

4.0 

0.5015 

1.6900 

0.0763 

0.2171 

0.0666 

1.1457 

5.0 

-0.0683 

-0.4060 

0.0042 

0.0120 

0.0017 

2.4701 


MMLE 

STABILITY 

& CONTROL 

DERIVATIVES 



CLA 

CMA 

CMQ 

CLDE 

CMDE 



2.7391 

-0.0887 

-0.1078 

0.5015 

-0.0683 




INITIAL 

INPUT DERIVATIVES 




2.4350 

-0.1089 

-0.6770 

1.6900 

-0.4060 
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B. pEst 

1. Model Validation 

A secondary purpose of this study was to integrate the 
pEst parameter estimation routine into the ongoing research 
program at the Naval Postgraduate School, Department of 
Aeronautics and Astronautics. Successful test cases for both 
the longitudinal and lateral-directional models were run to 
validate the model. The pEst parameter estimation routine 
should provide greater capability to the Department because it 
can be used to model the dynamic behavior of aircraft that 
cannot be appropriately modeled using the shorter PC-based 
routine. 

2. Longitudinal Model 

The longitudinal case was validated using a T-37 time 
history file received from NASA Dryden. The results of the 
validation are presented in Figures 16 through 19. The 
initial set up for the BQM-147 longitudinal model was 
presented in Table 9 in Chapter VI. The same current file was 
used for both the P2B1 and P2B2 maneuvers. 

a. P2B1 

The results of the P2B1 maneuver are shown in 
Figures 20 through 23 and the pEst output listing is shown in 
Table 12. The estimated response for a correlates well except 
at the extreme edge of the pitch-down part of the maneuver. 
The measured data shows a sinusoidal response at the beginning 
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of the window prior to any control input. This may be an 
accentuation of the process noise due to the low sample rate. 
On the initial run, the Cramer-Rao bound for C Lo was 
unacceptably high. Rather than accept the high bound, weight 
estimates based on discussions with NASA Dryden personnel, 
were applied to the outputs and subsequent runs showed 
improved results. The adjustment gave a better estimated 
response fit for a and a more reasonable Cramer-Rao bound, but 
the tradeoff was an increase in the pitch-damping coefficient 
which is evident in the q response. The response presented 
for both a and q represents the best possible fit achieved 
from multiple combinations of output weighting. The 
estimated response for normal acceleration correlates well, 
though it seems to lead the measured response by a time step. 
No time skewing was attempted. 
b. P2B2 

The results of the P2B2 maneuver are shown in 
Figures 24 through 27, and the pEst output listing is shown in 
Table 13. The final output weighting values used in the P2B1 
case were used for the P2B2 case because they appeared to give 
the best combination of results. The estimated a response 
does not keep up with the magnitude of the measured response 
in the pitch-down part of the maneuver and tends to dampen 
more quickly at the end of the maneuver. Again, this may be 
due to ignoring the 6 term in the a state equation. The 
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measured data show the same sinusoidal response in the absence 
of control input at the beginning of the time window and again 
at the end. The estimated q response correlates well. The 
damping on q tended to increase as weight on a was increased 
indicating a corellation between the a term and the q term. 
Normal acceleration fits well though it does not pick up the 
peaks at the edges of the maneuver. 

3. Lateral-Directional Model 

The lateral-directional case was validated using a PA- 
30 time history file received from NASA Dryden. The results 
of the validation are shown in Figures 28 through 33. The 
results of the combined P2B3-P2B5 lateral-directional maneuver 
are shown in Figures 34 through 39. The pEst output listing 
is presented in Table 14. The 8 response correlates well at 
the onset of the aileron input and follows the measured 
response through the rudder doublet. The r response fits well 
at the onset of the aileron input; however, the computed 
response appears to dampen more quickly than the measured 
response at the end of the rudder doublet. The p response 
does not follow the rate chrough the roll-left portion of the 
aileron doublet. The estimated response fits well during the 
rudder input but dampens more rapidly than the measured 
response. The roll angle response does not correlate well 
with the approximated angle data obtained through numerical 
integration of the roll rate. The lateral acceleration 
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response seems to match in phase but not in amplitude. The 
measured response shows better damping at the completion of 
the rudder doublet. 

4. Summary 

A summary of the longitudinal parameters is 
presented in Table 15. The columns labeled P2B1 and P2B2 
contain the results of the pEst parameter estimation runs. 
The Cramer-Rao bounds represent an approximation to the 
standard-deviation and are consistently one order of magnitude 
less than the value of the parameter. The large perturbation 
input maneuver and subsequent responses should lead one to 
suspect that the linear, PC based MMLE3 routine would be an 
insufficient model for this case study. The plots shown in 
Figures 12 through 15 support this suspicion. The pEst 
results provide a better fit between the measured and 
estimated data. The confidence level is higher for the 
parameters C L4e , Cm a and Cm 4e , which have lower relative Cramer- 
Rao bounds. The values for C La differ by a factor of 2 and 
the sign on Cm q was expected to be negative, indicating a 
restoring moment for positive pitch rate. The positive Cm q 
may indicate that there is a strong correlation between the a 
term and the q term. More data points are required to improve 
the confidence level for these two parameters. 

The lateral-directional case is summarized in Table 16. 
Running the lateral-directional case in the non-linear PC 
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based routine was unsuccessful because the process noise was 
too high. Only the pEst results are shown. The 6 derivatives 
and the p derivatives have reasonable Cramer-Rao bounds. The 
variance in the r derivatives is unacceptably high and 
indicates a low confidence level in the value of the 
parameters. A sign change in some of the control derivatives 
indicates a strong correlation with damping terms, a noted 
problem with parameter estimation methods. More flight data 
are required for any further analysis. 

Three options are available to correct for the high 
Cramer-Rao bound. 

1) Fix the value of the parameter 

2) Weight the estimate 

3) Use a better input stimuli 

Without a-priori knowledge of the parameters, fixing the 
value of the parameter was not a practical option. The 
weighting option did not help to improve the overall results 
because of the tradeoffs in the values of other parameters. 
A more reasonable approach to output weighting could have been 
made only if there were data other than the single flight test 
supporting the aerodynamic data base. With the limited 
maneuvers, better input stimuli were not available. Without 
any of the options available, the only choice is to accept the 
high Cramer-Rao bound with the understanding that the 
confidence level in the value of the parameter is very low. 
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There were three limitations in the data that affected the 
outcome of the lateral-directional results. Signal coverage 
is very important and unless the flight test maneuvers are 
restricted to small perturbations about a stable flight 
condition, the body angles should be included in the signal 
set. Noise in the same general freguency range as the 
response is a problem. The flight data were sampled at 16 Hz 
and nothing is known about the sensor accuracies, filtering, 
digitizing and time-tagging the data. The low sample rate may 
have distorted the response signal by accentuating the process 
noise. To excite the lateral-directional modes, rudder and 
aileron doublets should be included in a single input 
maneuver, separated by sufficient time to allow for the first 
response to dampen. In Flight 03-24, the maneuvers were 
independent and combining the data sets together may have 
induced additional error in the response estimates. 
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Figure 17 pEst Validation: a Response 
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Figure 21 BQM-147 a Response (P2F1) 
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Figure 23 BQM-147 Normal Acceleration (P2B1) 
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TABLE 12 pEst OUTPUT LISTING 


Iteration 

tev-Marq 

17 

used 1 trials. 

v • 0. 

len » 

3.150 

coat « 37.65 

name 

value on? 

delta 

name 

value 

on? 

delta 

cNormO 

1.173 T - 

•4573E-02 

cmq 

0.3234 

T 

0.2962E-02 

cNorna 

0.8794E-01 T 0 

. 1 699E-03 

cmde 

-.1466E- 

•02 T 

0.4741E-05 

cNormde 

0.2389E-01 T - 

.3282E-03 

alphaO 

-.5625 

T 

-. 1678E-02 

cmO 

-.3865E-01 T - 

-1539E-05 

ka 

1.663 

T 

0.8075E-02 

cma 

-.2951E-02 T - 

•1587E-04 






total cost « 37.65 

cost per response : 

alpha q an 

35.5 69.5 7.88 

Iteration 17 used 12 Integrations 
••• iteration converged 
estimate used 209 Integrations 


pEst: par 
name 

♦cr 

value on? 

cr bound 

name 

value 

on? 

cr bound 

cNormO 

1.173 

T 

0.4632E-01 

cmq 

0.3234 

T 

0.2749E-01 

cNorma 

0.8794E-01 

T 

0.3379E-02 

cmde 

- . 1466E- 

-02 T 

0.4138E-04 

cNormde 

0.2389E-01 

T 

0.2481E-02 

alphaO 

- .S625 

T 

0.1954 

cmO 

-. 3865E-01 

T 

0.1263E-02 

ka 

1.663 

T 

0.6568E-01 

cma 

pEst: 

- .2951E-02 

T 

0.1338E-03 
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Figure 25 BQM-147 a Response (P2B2) 
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Figure 27 BQH-147 Normal Acceleration (P2B2) 
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TABLE 13 pEst OUTPUT LISTING (P2B2) 


Iteration 14 


Lev-Mare a 

sed 1 trials 

. v « 0. 

len * 

1.637 cost • 103.0 

name 

value on? 

delta 

name 

value on? 

delta 

cNormO 

0.6941 T 

-.1956E-02 

cma 

-.2266E-02 T 

- .3263E-05 

cNorma 

0.4120E-01 T 

-.1341E-03 

croq 

0.1355 T 

0.3866E-02 

cNormde 

0.2149E-01 T 

-.5021E-04 

cmde 

-.1600E-02 T 

0.1210E-04 

cmO 

-.3293E-01 T 

0.2963E-05 

ka 

1.573 T 

- .4011E-04 

total cost 

« 103.0 






cost per response : 

alpha q an 

66.9 223. 19.5 

Iteration 14 used 11 Integrations 
••• iteration converged 
estimate used 157 Integrations 

pEst: write out2.p2b2 

writing computed time history to file *out2.p2b2* 
pEst: par +cr 


name 

value on? 

cr bound 

na*e 

value on? 

cr bound 

cNormO 

0.6941 

T 

0.2587E-01 

CM 

-. 2266E-02 

T 

0.8947E-04 

cNorma 

0.4120E-01 

T 

0.1786E-02 

cwq 

0.1355 

T 

0.1741E-01 

cNormde 

0.2149E-01 

T 

0.1869E-02 

erode 

-.1600E-02 

T 

0.4203E-04 

cmO 

--3293E-01 

T 

0.1119E-02 

ka 

1.573 

T 

0 3094E-01 
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Figure 31 pEst Validation: p Response 
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Figure 34 BQM-147 Lateral-Directional Input 
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Figure 35 BQM-147 8 Response 
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Figure 37 BQM-147 p Response 
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Figure 38 BQM-147 <p Response 
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TABLE 14 pEst OUTPUT: LATERAL-DIRECTIONAL 


Lev-Marq 

used 1 trials, v » 0. 

len - 

3.436 cost * 1003. 

name 

value on? 

delta 

raise 

value on? delta 

cyO 

0.4222E-02 

T 

- .5625E-06 

cldr 

0.38S6E-04 

T -.2546E-06 

cyb 

-. 8251E-02 

T 

-.1725E-04 

cnO 

0.3708E-03 

T -2122E-05 

cydr 

-. 1862E-02 

T 

0.2043E-05 

cnb 

0.4335E-03 

T 0.3397E-06 

clO 

0.1080E-02 

T 

0.3836E-0S 

cnp 

- .1212E-01 

T 0.1558E-03 

clb 

-.9698E-04 

T 

0.2407E-07 

cnr 

-.7893E-03 

T 0.1786E-03 

clp 

-. 5860E-01 

T 

- .2733E-03 

cnda 

-.7178E-04 

T 0.3820E-06 

clr 

-. 3279E-02 

T 

- .1693E-03 

crdr 

-.2109E-03 

T 0.2517E-06 

clda 

- . 6099E-03 

T 

-.1818E-0S 

kb 

1.067 

T 0.2492E-02 


total coat - 1003. 

cost per response : 

beta p r phi 

51.2 3.714E*03 322. 817 

iteration 8 used 19 integrations 
* * * Iteration converged 
estimate used 155 Integrations 

pEet: par *cr 


name 

value on? 

cr bound 

name 

value on? 

cr bound 

cyO 

0.4222E-02 

T 

0.272SE-02 

cldr 

03856E-04 

T 

0.8644E-05 

Cyb 

- . 3251E-02 

T 

0.1162E-02 

cnO 

0.3708E-03 

T 

0.9765E-04 

Cydr 

- . 1882E-02 

T 

0.6809E-03 

cnb 

0.4335E-03 

T 

0.1105E-04 

CIO 

0.1080E-02 

T 

0.4500E-04 

cnp 

- .1212E-01 

T 

0.3082E-02 

Clb 

- 9698E-04 

T 

0.9580E-0S 

cnr 

- .7893E-03 

T 

0.6030E-02 

clp 

- .5860E-01 

T 

0.2814E-02 

cnda 

- .7178E-04 

T 

0.1808E-04 

clr 

--3279E-02 

T 

0.2653E-02 

cndr 

- .2109E-03 

T 

0.1095E-04 

clda 

-.6099E-03 

T 

0.203 3E-04 

kb 

1.067 

T 

0.9169E-01 


pEst: write outlat7 


ay 

109. 
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TABLE 15 LONGITUDINAL SUMMARY 


Parameter 

Initial 

MMLE3 

P2B1 

P2B2 

c La 

2.4350 

2.7391 

5.0367 

2.3600 

CL* e 

1.6900 

0.5015 

1.3689 

1.2314 

Cm a 

-0.1089 

-0.0887 

-0.1691 

-0.1298 

Cm q 

-0.6670 

-0.1078 

0.4997 

0.1355 


-0.4060 

-0.0683 

-0.0840 

-0.0917 1 

Cramer-Rao Bounds 

c La 


0.1084 

0.1936 

0.1023 



0.0763 

0.1422 

0.1071 

cm a 


0.0028 

0.0077 

0.0051 

cm q 


0.0414 

0.0275 

0.0174 



0.0042 

0.0024 

0.0024 
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TABLE 16 LATERAL DIRECTIONAL SUMMARY 


Parameter 

Initial 

pEst 

Cramer-Rao 

C V8 

-0.2005 

-0.4728 

0.0665 

c i0 

-.0325 

-0.0056 

0.0005 

Cn fi 

0.0430 

0.0248 

0.0006 

c1 p 

-0.3525 

-0.0586 

0.0028 

Cn p 

0.0027 

-0.0121 

0.0031 

Cl r 

0.0729 

0.0033 

0.0026 

Cn r 

-0.0230 

-0.0008 

0.0060 

ci 5a 

0.1377 

-0.0349 

0.0012 

Cn 5a 

-0.0194 

-0.0041 

0.0010 

C Yir 

0.1130 

-0.1078 

0.0390 

Cl« r 

0.0023 

0.0022 

C.0005 

Cn «r 

-0.0138 

-0.0121 

0.0006 



















VIII. CONCLUSIONS AND RECOMMENDATIONS 


A. CONCLUSIONS 

Based on the results presented in Chapter 7, the following 
conclusions are made. 

1) The PC-based parameter estimation routine is ideally 
suited for small maneuvers on stable platforms. The 
magnitude of the maneuvers in Flight 03-24 exceeded the 
range of validity for the linear model. 

2) The pEst routine was successfully integrated into the 
NPS Department of Aeronautics and Astronautics flight 
research program. 

3) Noise in the same general frequency range as the system 
response constitutes a major problem. 

4) Body angles should be included in the signal set unless 
the maneuvers are restricted to small perturbations about a 
stable flight condition. 

5) Improving the confidence level in the values of the 
parameters presented in Tables 15 and 16 will require 
modifications to the flight test program. 


B. RECOMMENDATIONS 

Based on the conclusions presented above and the problems 
encountered during the course of this study, the following 
recommendations are made. 

For the Marine Corps UAV Project Office, Intelligence 
Systems (C2IU), Marine Corps Research, Development and 
Acquisition Command, Quantico, VA. : 
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1) Evaluate the current status of the BQM-147 Program and 
determine if additional flight testing to provide higher 
confidence values of stability and control parameters is 
warranted. 

2) If higher confidence values are desired, formulate a 
clear set of objectives and modify the flight test program 
to include a series of multiple maneuvers designed 
specifically to meet those objectives. The maneuvers should 
include repeated data points, and roll and pitch angles 
should be measured. 


For the NPS, Department of Aeronautics and Astronautics: 

1) Update the Department PC Lab to include MATLAB version 
3.5j or higher. As a part of the update, add the State- 
Space Identification Tool to the current inventory of MATLAB 
toolboxes. The update will allow instructional use of the 
PC-based routine in the flight dynamics and systems courses 
and practical application of a desk-top routine at the NPS 
UAV Flight Research Lab. 

2) Incorporate pEst as an instructional tool in the Flight 
Dynamics courses at NPS. 

3) Investigate modification of the user-modifiable 
subroutines in pEst to include helicopter dynamics and 
missile dynamics. 

4) Investigate the feasibility of offering a course in non¬ 
linear flight dynamics or non-linear control of aerospace 
vehicles. 

5) Suggest as a potential project for the Avionics System 
Design or Aircraft Design courses, the design of an 
instrumentation package for a UAV Flight Test Program that 
would respond to a real-world DoD requirement. 
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Figure A-l EXDRONE Baseline Configuration 
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APPENDIX B FLIGHT 03-24 TIME HISTORY DATA 


A. PITCH-UP DOUBLET (P2B1) 



Figure B-2 Angle of Attack Response 
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PITCH-DOWN DOUBLET 



re B-5 Elevator Control Input 



Figure B-6 Angle of Attack Response 
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C. COMBINED ROLL AMD YAW DOUBLET (P2B3-P2B5) 




Figure B-10 Sideslip Response 
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Lateral Acceleration (g) 



Figure B-13 Lateral Acceleration 
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